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Abstract In this contribution we discuss about numerical modeling of nebulae. In 
particular we emphasize on the dynamical evolution of an H II region and the chem- 
ical structure of a Photodissociation region. We do this using the Smoothed Parti- 
cles Hydrodynamics code SEREN and the recently developed astrochemistry code 
^ (-h \ 3D-PDR, respectively. We show an example application by simulating a cometary 

O I globule using these two codes. 



1 Introduction 

Nebulae (H II regions) are large regions consisting of ionised gas, particularly of 
hydrogen. The source of ionisation is usually a single or multiple massive stars 
emitting ultraviolet radiation with photons carrying more energy than the ionisation 
potential (hv > 13.6eV). A nebula is mainly structured by three different parts; the 
ionised region (number density of n < 200cm -3 and temperature of T ~ 10 4 K), 
the photodissociation regiorQ (PDR; 200 < n < 10 5 cm~ 3 , 10 - 20 < T < 10 4 K), 
and the dark molecular region (n> 10 5 cm -3 , T ~ 10 — 20 K). Of particular interest 
are PDRs; they are ubiquitously present in the interstellar medium (ISM) consisting 
of predominantly neutral gas and dust illuminated by FUV radiation (6 < hv < 
13.6eV) and they occur in any region of the ISM that is dense and cold enough to 
remain neutral but has too low column density to prevent the penetration of FUV 
photons. Over the past few decades, effort has been made to study numerically the 
physics and chemistry of nebulae. Due to computational speed and memory capacity 
issues, detailed three-dimensional simulations are a reality since the past 5-7 years 
only. Even nowadays the numerical codes are divided in two main categories; codes 
that study the dynamical evolution and codes that study the chemical structure. An 
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integrated code that treats simultaneously detailed dynamics, UV propagation and 
chemistry offering a realistic temperature and therefore pressure structure in H II 
regions is still lacking (although significant effort in this direction has been made by 
E] E] El). 



2 Dynamical & Chemical modeling 

Several techniques have been proposed by various workers on dynamical modeling 
of H II regions in adaptive mesh refinement codes and in smoothed particle hydro- 
dynamics codes. have proposed a HE ALP ix-based @ algorithm to simulate the 
propagation of the UV radiation in the ISM by invoking the on-the-spot approxi- 
mation |[T3l . This algorithm has been incorporated in the SPH code SEREN and it 
creates hierarchies of rays emanated spherically symmetric from the excited source. 
Each ray splits up to four child-rays wherever the resolution of the radiation trans- 
fer matches the resolution of the hydrodynamics, locally. This algorithm adopts two 
equations of state; an isothermal for the ionised region (T = 10 4 K everywhere), and 
a barotropic for the neutral region (see Eqn. 6 [3 ]). The temperature between these 
two regions, and therefore where the PDR is located, is smoothed using a linear in- 
terpolation. Although this technique speeds up the thermodynamical calculations, it 
significantly limits the ability to study the chemistry of PDRs which is quite impor- 
tant to understand in order to explore the insights of the chemical structure and star 
formation in ionised nebulae. We thus need further tools to implement. 

Perhaps the most challenging part of modeling a nebula is its chemistry. That 
is because one has to take into account a realistic three-dimensional treatment of 
the UV radiation (i.e. as obtained by MOCASSIN EEl), self- shielding of individ- 
ual species against the UV radiation, cooling and heating processes and a com- 
plicated network of reactions. Of particular interest is the three-dimensional mod- 
eling of PDRs which has been achieved in 3D -PDR for the first time. 3D -PDR 
solves the chemistry and the thermal balance self-consistently within a given three- 
dimensional cloud of arbitrary density distribution. It uses the chemical model 
features of the fully benchmarked one-dimensional code UCL_PDR [lj and a ray- 
tracing scheme based on the HEALPix package to calculate the total column den- 
sities and thus to evaluate the attenuation of the FUV radiation into the region, and 
the propagation of the FIR/submm line emission out of the region. An iterative cycle 
is used to calculate the cooling rates using a three-dimensional escape probability 
method, and heating rates. At each element within the cloud, it performs a depth- 
and time- dependent calculation of the abundances for a given chemical network 
to obtain the column densities associated with each individual species. The itera- 
tion cycle terminates when the PDR has obtained thermodynamical equilibrium, in 
which the thermal balance criterion is satisfied i.e. the heating and cooling rates are 
equal to within a user-defined tolerance parameter. 3D -PDR determines the relative 
abundances of a limited number of atomic and molecular species at each cloud ele- 
ment, by solving the time-dependent chemistry of a self-contained network of for- 
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Fig. 1 Top row from left to right: emission maps of [C II] 158 fim, [C I] 610 jim and [O I] 63 jim. 
Bottom row from left to right: emission map of CO (1-0), RGB composite image with colour bar 
ratios of 8:1:2 for CO(1-0):[C i]:[C II], cross-section density profile at z = Opc. 



mation and destruction reactions. We use the UMIST database containing 33 species 
(including e~) and 320 reactions. We solve for steady-state chemistry (chemical 
evolution time set to t = lOOMyr), although the code is able to follow the full time 
dependent evolution within the cloud. 

This treatment of PDR chemistry offers a more accurate temperature profile in 
comparison with the approximations made in the dynamical simulations. However, 
due to the very high computational demanding, it is impossible to include such 
complicated UV and PDR calculations in a dynamical code, unless new techniques 
are implemented. 



3 Example: simulating a cometary globule 

In this example we perform two different simulations using the codes described 
above. In the first simulation (SPH run) we use the SEREN code to evolve dy- 
namically an initially uniform density spherical clump (radius R = 0.5 pc, mass 
M = 20 M , temperature T = 100K, consisting purely of atomic hydrogen) as it 
interacts with an external and approximately plane-parallel radiation field. In the 
second simulation (PDR run) we use the 3D -PDR code in one of the snapshots 
selected from the SPH run in which we examine the PDR chemical structure. 

In the SPH run, as the clump interacts with the external radiation field emitted 
spherically by an excited source (at a rate of ^{ yC = 3.2 x 10 48 photons s _1 , placed 
D = 3.5 pc away from the centre of the clump) the ionising radiation boils off its 
outer layers from the side in which the flux is impinging. The shock front that is 
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formed compresses the remaining neutral gas turning it into a rod shape. However, 
the internal thermal pressure of the neutral gas is constantly increasing and this 
results in its re-expansion (see §4.4 of O for full discussion), which may lead to 
star formation |[T0l [3ll. 

In the PDR run, we take a snapshot from the SPH run at t = 0.12Myr and we 
use it as initial conditions in 3D -PDR. A cross section plot of the density structure 
at that time is shown at the bottom right of Fig. [T] where the ionising radiation is 
impinging from bottom to top. The UV photon flux corresponds to a field strength 
of approximately % = 50Draines. The top three and the bottom left frames show 
the emission maps for the most dominant coolants. The middle frame at the bottom 
shows an RGB composite image of three different emission maps. From these maps 
we see that the species in the PDR are distributed smoothly and follow the density 
profile. The weakest emission is produced by [O I] 63 jim and the strongest by 
CO (1-0) implying that the molecular gas dominates over the atomic contribution. 
However, considering the transition frequencies for these maps, we find that the 
[O I] 63 fim and the [C II] 158 fim lines are the dominant coolants. Full details and 
results are discussed in H. 



4 Conclusions 

We discussed about the dynamical and chemical modeling of H II regions using the 
SEREN SPH code and the recently developed 3D -PDR code respectively. We per- 
form an example of a cometary globule and we show results from these two codes. 
Since we lack of algorithms treating detailed dynamical and chemical calculations 
in ionized regions, effort has to be made towards this direction. 
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